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We investigate the clustering and phase separation of a model of ultrasoft, oppositely charged 
macroions by a combination of Monte Carlo and Molecular Dynamics simulations. Static and dy- 
namic diagnostics, including the dielectric permittivity and the electric conductivity of the model, 
show that ion pairing induces a sharp conductor-insulator transition at low temperatures and densi- 
ties, which impacts the separation into dilute and concentrated phases below a critical temperature. 
Preliminary evidence is presented for a possible tricritical nature of the phase diagram of the model. 
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Clustering of oppositely charged ions is a common 
mechanism [TJ [2] which strongly influences the thermo- 
dynamic and transport properties of electrolytes, in par- 
ticular the phase separation into dilute and concentrated 
ionic solutions predicted to occur at low temperatures [3] . 
This phase separation (akin to vapour-liquid coexistence) 
has been investigated in considerable detail both theoret- 
ically and by computer simulations in the context of the 
primitive model of electrolytes, consisting of oppositely 
charged hard spheres immersed in a dielectric continuum 
(implicit solvent), with most of the published work fo- 
cusing on the "restricted" version of the model (RPM) 
featuring equal anion and cation diameters [3H7]. The 
RPM is a reasonable model for strong, monovalent elec- 
trolytes, like aqueous solutions of NaCl. The formation 
of long-lived dipolar pairs considerably affects the coex- 
istence curve Hj , and severely limits the ergodicity of 
traditional Monte Carlo (MC) or Molecular Dynamics 
(MD) simulations at low temperatures [5]. 

Recently the attention has shifted to the rich phase be- 
haviour of "colloidal electrolytes" where the anions and 
cations are highly charged, hard colloidal [TU] or nano- 
metric particles [IT], usually in the presence of added 
salt H2HH]. In this Letter, we extend the RPM to a 
broader class of soft matter systems by introducing an 
"ultrasoft" restricted primitive model (URPM), where 
macroions are assumed to be penetrable charged par- 
ticles, i.e., the hard cores of charged colloids are re- 
placed by bounded interactions at short interionic dis- 
tances. Ultrasoft core representations of the effective 
interaction between the centres of mass (CM) of poly- 
mer coils have proved very successful in describing dilute 
and semi-dilute polymer solutions [T5J [T5]. Our model 
generalises such a representation to solutions of oppo- 
sitely charged polyelectrolytes chains [TTJ [TB]. We ex- 
plore, through MC and MD simulations, the subtle in- 
terplay between the clustering effects associated with in- 
terpenetrating soft core particles [19] and the long-range 
Coulombic interactions, and its influence on phase sep- 



aration. Our simulation results for the URPM reveal a 
non-trivial topology of the phase diagram of the model 
and suggest a strong link between phase separation and 
a purely classical conductor-insulator (CI) transition at 
low temperatures, reminiscent of the behaviour of liquid 
metals [2U]. 

The URPM is a system of n+ cations of total charge 
+Q and n_ anions of charge — Q per unit volume moving 
in a dielectric continuum of relative permittivity e'; global 
charge neutrality implies n + = n_ = n. The charge 
around the CM of each macroion follows a Gaussian dis- 
tribution Q a pa(s) — Q a exp[— S 2 /2(7 2 ] (27TCT 2 ) 3 ^ 2 (a — 

+ or -) where s is the distance from the CM. The result- 
ing pair potentials as functions of the distance r between 
the CM's of two ions are: 
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where — uq = — Q 2 / ^/ne'a is the energy of a fully overlap- 
ping anion/cation pair (r — 0). For r > a, the pair poten- 
tials go over to the Coulombic interaction between point 
charges. A related model was used previously to inves- 
tigated a very different system, namely a semi-classical 
Hydrogen plasma under astrophysical conditions of high 
temperatures and densities 21j. In the following, we 
will use y/2a, Uq, and yjlma 2 /u as units of length, 
energy and time, respectively. We will report MC and 
MD simulations results over a broad range of total den- 
sity p = 2n and temperature T, obtained for samples of 
N = N + + iV_ = 1000 anions and cations under periodic 
boundary conditions, employing the Ewald summation 
technique. Comparison of different simulations methods 
and different thermal histories provide evidence of proper 
equilibration of our samples over the investigated range 
of state parameters. 

The classical ground state energy (T = 0) is U = 
—Nu , which is extensive, ensuring the existence of a 
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FIG. 1: Radial distribution functions for selected tempera- 
tures (as indicated) along the isochore p = 0.01: (a) (r), 

(b) 9++(r). 



proper thermodynamic limit |22) . In the high tempera- 
ture limit (T > 0.5) we found that the model is accurately 
described by the random phase approximation (RPA), 
which reduces to the familiar Debye-Hiickel theory for 
point ions and has been shown to be very accurate also 
for other soft core models [TS1 [THl dH] at high densities. 

An interesting phenomenology appears at sufficiently 
low densities and temperatures, where the system shows 
clear signatures of clustering. To quantify this phe- 
nomenon we first inspect the equilibrium pair structure, 
which is characterised by the two radial distribution func- 
tion (RDF) g ++ (r) = g (r) and g+-(r). MD data 

are shown in Fig. [T] along the isochore p = 0.01, rep- 
resentative of the behaviour of the system at low den- 
sity. At low T the sharp rise of g^ (r) as r — >• points 

to strong anion/cation clustering, as evident in snap- 
shots of ion configurations. Furthermore, clustering at 
low T is also indirectly evident from g ++ (r), which ap- 
pears to imply an attraction between equally charged ions 
[9++( r <t ) > 1] — a f ac t that can only be rationalised 
by the formation of tight neutral pairs or higher order 
clusters. 

From the Fourier transforms of the pair correlation 
functions h a p(r) = g a /3(r) — 1 (a, j3 = + or -), we ex- 
tract the charge-charge structure factor 



Scc(fc) = l + n[h ++ (k)-h+_(k)} 
k 2 



k 2 + exp[— k 2 a 2 ] 
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where kd = (87mQ 2 / V k^T^j is the inverse Debye 
screening length, k is the wave-number and hats im- 
ply Fourier transforms. The RPA incorporates the exact 
Stillinger-Lovett sum rule [2] lim^o K^Scc(k)/k 2 = 1, 
valid for a conducting medium. The sum rule is well 
satisfied by our simulation data for Scc{k) along the iso- 
chore p = 1 at all temperatures, and the RPA expres- 
sion [Eq. (12])] provides an accurate representation of the 
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FIG. 2: Charge-charge structure factors Scc(k) for selected 
temperatures (as indicated) along the isochore p = 0.01 from 
MD simulations (filled symbols) and from Eq. H (solid lines). 
The dashed line indicates a quadratic fit (fc/K) 2 , where k is 
an effective inverse screening length, to the low-fc behaviour 
of Sec (k) for T = 0.025. 



data for all k at sufficiently high T (not shown here). 
Along the low density isochore p — 0.01 (see Fig. [21), 
the RPA expression is accurate only for T > 0.5 while 
for T < 0.1, strong deviations from the perfect screening 
sum rule are observed; fits to the small-fc parabolic be- 
haviour of Scc{k) lead to an effective inverse screening 
length k systematically larger than acd (by a factor of two 
at T = 0.025). This implies that the URPM no longer 
behaves as a conductor but rather as a dielectric medium 
of neutral anion/cation clusters. 

Using a standard geometric definition of n-ion clusters, 
namely that n ions form an n-mer if each ion lies within 
a distance r c (= 1.0) of at least one other ion in the clus- 
ter, we have determined the percentages P n of n-mers, 
averaged over all configurations. Examples of Pi (iso- 
lated ions), P 2 (anion/cation pairs) and P4 (tetramcrs) 
as functions of T are shown in Fig.[3]-a, along the isochore 
p = 0.01. As expected, P2 is close to 100% at the lowest 
T, and drops rapidly for T > 0.03, while Pi starts from 
and increases towards 100% for T > 0.03. Note that 
the percentage of trimers (now shown) is always negligi- 
ble while the percentage of tetramers can be significant 
(> 5% at the lowest temperatures). The lifetime r 2 of 
pairs, as estimated from MD simulations, increases dra- 
matically as T drops below 0.05, and is approximately 
two orders of magnitude larger than that of the other 
n-mers. 

To assess explicitly the effect of pairing on the dielec- 
tric properties of the system, we calculate the dielectric 
permittivity e from the fluctuations of the total dipolc 
moment of the simulated periodic sample, M — QiYi, 
according to the standard Kirkwood relation [23] , Re- 
sults from our simulations along the isochore p = 0.01 
are shown in Fig. [3]-b. For T < 0.03 the fluctuations of 
M are small, and the resulting e takes on values of the 
order of 1.3, typical of dielectrics made up of polarizable 
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FIG. 3: Comparison of (a) percentages of selected n-mers 
P n , (b) dielectric permittivity e, and (c) conductivity <t c as 
functions of T along the isochore p = 0.01. 



molecules (ion pairs in our case). At higher T fluctua- 
tions of M are strongly enhanced, due to the break-up of 
ion pairs, and e rises sharply towards values larger than 
100, typical of a conductor (for an infinite conducting 
sample, e — > oo). 

These suggestions of a CI transition prompted us 
to extract some dynamic diagnostics from MD simula- 
tions [29] • Figure [4] shows plots of the time-dependent 
electric dipole diffusion c(t) = (|M(t) - M(0)| 2 ). For a 
conductor, c(t) increases linearly at long times (accord- 
ing to the generalized Einstein relation) , and the asymp- 
totic slope is proportional to the electrical conductivity 
a c . At the lowest T, the slope vanishes, i.e., a c = 0, 
corresponding to a dielectric insulator state. The re- 
duced angular frequency of the oscillations observed at 
low T is close to the harmonic oscillator frequency of the 
parabolic anion/cation potential well [Eq. Q], confirm- 
ing the formation of long-lived pairs that do not con- 
tribute to the electrical conductivity. The temperature 
variation of a c is illustrated in the inset of Fig. [4] along 
two isochores and confirms the sharp drop of er c by sev- 
eral orders of magnitude over a narrow range of temper- 
atures. The conductivity data have been fitted by power 
laws A(T ~ Tu) 1 - 2 , yielding the apparent CI transition 
temperatures T a — T a (p). 

Figure [3] provides an overview of the behaviour of the 
system along the representative isochore p = 0.01, dis- 
playing results for the percentages P n of n-mers, the per- 
mittivity e, and the conductivity cr c as functions of T. 
The three plots show a strong correlation between the 
sharp rise in e, a c and Pi around T = 0.03, indicative of 
a CI transition, driven by pairing and rounded by finite 
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FIG. 4: Diffusion of the total dipole moment c(t) = (\M(t) - 
M(0)| 2 ) as a function of reduced time t for selected tempera- 
tures (as indicated) along the isochore p = 0.01. Inset: elec- 
trical conductivity cr c as a function of T for p = 0.006 (filled 
circles) and p — 0.02 (empty squares). Solid lines are power 
law fits (7 C (T) w A(T - T a ) 1,2 . 



size effects. 

The final step of our investigation is to relate the clus- 
tering to the liquid-vapour phase transition expected at 
low T. To that purpose we have carried out grand- 
canonical MC simulations in periodic cells of side L = 
9.28 and L — 14.7, combining biased MC sampling and 
histogram reweighting techniques [24] to construct the 
coexistence curves shown in Fig. [5] together with the 
loci of state points where the CI transition is expected 
(i.e., a c — > according to the aforementioned power law 
fits) and where the fraction Pi of free ions reaches 30%. 
Different choices of the value of Pi at the pairing tran- 
sition temperature, as well as variations of r c by some 
10% or 20%, give rise to pairing transition temperatures 
shifted by at most ±0.01 with respect to those displayed 
in Fig. [5] The extrapolation of the two transition lines 
intersect the coexistence curve close to the critical point, 
roughly estimated to be T c ~ 0.018, p c ~ 0.05 on the ba- 
sis of the results for L — 14.7. The phase diagram shows 
that the liquid-vapour coexistence curve is better fitted 
by a scaling law with a critical exponent /? = 1 than 
by the Ising universality class exponent ft = 0.326 (or 
by the mean- field exponent ft = 0.5, not shown). Inclu- 
sion of corrections to scaling did not improve the fit for 
ft = 0.326 substantially. The break-down of an Ising-like 
description of the coexistence curve, the apparent value 
of the fitted critical exponent (ft « 1) and the proximity 
of the critical point to the CI line hints at a tricritical 
nature of the phase diagram of the URPM, very different 
from the Ising-like phase diagram of the RPM. This con- 
jecture will be tested by finite size scaling calculations 



4 



0.04- 



0.03 



0.02 



0.01 



«4 



Or 



e • 



c » 



b O 
en • 



Coexistence (L=9.3) 
Coexistence (L=14.7) 
Pairing 
P=1.0 
B=0.326 




0.01 



0.1 



FIG. 5: Extended phase diagram of the URPM: liquid-vapor 
coexistence from GCMC simulations (L — 14.7, filled circles; 
L = 9.3, empty circles), pairing transition points (T p ,p p ) de- 
termined from the condition P\{T p ,p p ) — 30% (empty trian- 
gles) and CI transition points (T CT ,p CT ) (stars; see text for def- 
inition). Critical fits to the coexistence curve for L = 14.7 are 
shown for /3 = 0.326 (dotted line) or /3 = 1.0 (solid line). The 
estimated critical point corresponding to f3 = 1.0 is shown as 
a bigger filled sphere. The law of rectlinear diameter is also in- 
cluded (dash-dotted line and filled squares). The dashed line 
through the pairing transition points is drawn as a guide to 
the eye. Inset: MD-generated configuration at p = 0.01 and 
T = 0.00125. White and red spheres (not drawn to scale) 
indicate oppositely charged particles. 



tween the hard core and soft core models may be traced 
back to the fundamental difference between long-lived 
ion pairs at low temperatures and densities: these pairs 
are strongly dipolar dumbbells in the RPM, while in the 
present URPM they are non-polar, polarizable entities. 

In summary, we have introduced a simple model of 
oppositely charged, interpenetrating macroions, which 
generalizes the familiar RPM to "soft" polyelectrolytes. 
Using static and dynamic diagnostics in MC oand MD 
simulations, we have provided a quantitative character- 
ization of pairing and clustering of ions at low temper- 
atures and densities, and their impact on the segrega- 
tion into coexisting dilute and concentrated phases. The 
predicted clustering and segregation of the URPM are 
reminiscent of the experimentally observed complexation 
of anionic and cationic polyelectrolytes and subsequent 
complex coacervation |17| 118] , as confirmed by approxi- 
mate field-theoretic calculations and simulations [37J [55] . 
Examination of the state-dependence of the dielectric 
pemittivity and of the electric conductivity suggests that 
pairing leads to a CI transition resembling that observed 
in liquid metals, such as Hg or Rb. Future work will 
concentrate on a quantitative analysis of finite size ef- 
fects and the extension of the URPM to the unrestricted, 
asymmetric case. 

We thank C. Valeriani, K. Binder, and S. Clarke for 
useful discussions. D.C. and G.K. acknowledge financial 
support by the Austrian Science Foundation (FWF) un- 
der Proj. No. P19890-N16. 



in future work. We note that the putative tricritical be- 
havior of our model may be different in nature from the 
one observed in a lattice version of the RPM (see [S3] 
and references therein). The latter behavior is linked in 
fact to an order-disorder transition, reminiscent of that 
observed in some antiferromagnetic materials. 

As in the case of the RPM, the dilute phase is dom- 
inated by long-lived ion pairs (and larger neutral clus- 
ters), while in the dense phase, ions are essentially free. 
This is illustrated by the snapshot in the inset of Fig. [5] 
which shows droplet formation during an MD simulation 
at low density and temperature. Identifying the length 
scales a of the RPM and the URPM, the reduced crit- 
ical parameters T c and p c of the two models are of the 
same order of magnitude, suggesting that the phase sep- 
aration is essentially of Coulombic origin. However, a 
closer comparison between Fig. [5] and recent simulation 
results for the RPM (25] reveals a striking difference in 
the topology of the phase diagram: while in our model 
the pairing transition line intersects the coexistence curve 
close to the estimated critical point, pair formation in the 
RPM represents a crossover occurring at densities much 
lower than the critical one. Thus the discrepancies be- 
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